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Abstract 

We review and apply Quasi Monte Carlo (QMC) and Global Sensitivity Analysis (GSA) tech¬ 
niques to pricing and risk management (greeks) of representative financial instruments of in¬ 
creasing complexity. We compare QMC vs standard Monte Carlo (MC) results in great detail, 
using high-dimensional Sobol’ low discrepancy sequences, different discretization methods, 
and specific analyses of convergence, performance, speed up, stability, and error optimisation 
for finite differences greeks. We find that our QMC outperforms MC in most cases, including 
the highest-dimensional simulations and greeks calculations, showing faster and more stable 
convergence to exact or almost exact results. Using GSA, we are able to fully explain our 
findings in terms of reduced effective dimension of our QMC simulation, allowed in most 
cases, but not always, by Brownian bridge discretization. We conclude that, beyond pricing, 
QMC is a very promising technique also for computing risk figures, greeks in particular, as 
it allows to reduce the computational effort of high-dimensional Monte Carlo simulations 
typical of modern risk management. 
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1 Introduction 


Nowadays market and counterparty risk measures, based on multi-dimensional, multi-step Monte 
Carlo simulations, are very important tools for managing risk, both on the front office side, for 
sensitivities (greeks) and credit, funding, capital valuation adjustments (CVA, FVA, KVA, gener- 
ically called XVAs) and on the risk management side, for risk measures and capital allocation. 
Furthermore, they are typically required for regulatory risk internal models and validated by 
regulators. The daily production of prices and risk measures for large portfolios with multiple 
counterparties is a computationally intensive task, which requires a complex framework and an 
industrial approach. It is a typical high budget, high effort project in banks. 

In the past decades, much effort was devoted to the application of Monte Carlo technique^] 
to derivatives pricing |Boy77j IBBG971I.TacOll iGl af)3] . The main reason is that complex financial 
instruments usually cannot be priced through analytical formulas, and the computation of high¬ 
dimensional integrals is required. Monte Carlo simulation is, then, a common way to tackle 
such problems, since it reduces integration to function evaluations at many random points and 
to averaging on such values. As a result, virtually any product can be easily priced in any 
dimension. However, this method is rather time consuming and the convergence rate is slow, 
since the root mean square error decays as V -1 / 2 , where N is the number of sampled points. 
Various “variance reduction” techniques exist, which can improve the efficiency of the simulation, 
but they don’t modify the convergence rate [■lacOll iGlaM] , 

Quasi Monte Carlo represents a very efficient alternative to standard Monte Carlo, capable 
to achieve, in many cases, a faster convergence rate and, hence, higher accuracy |.Tac011 [GlaMl 
IMF991 ISKOSbl ISKOhal IWa.nfffll IKFSM111ISAKK12] . The idea behind Quasi Monte Carlo meth¬ 


ods is to use, instead of pseudo-random numbers (PRN), low discrepancy sequences (LDS, also 
known as quasi-random numbers) for sampling points. Such LDS are designed in such a way 
that the integration domain is covered as uniformly as possible, while PRN are known to form 
clusters of points and always leave some empty areas. Indeed, the very random nature of PRN 
generators implies that there is a chance that newly added points end up near to previously 
sampled ones, thus they are wasted in already probed regions which results in rather low con¬ 
vergence. On the contrary, LDS “know” about the positions of previously sampled points and 
fill the gaps between them. Among several known LDS, Sobol’ sequences have been proven 
to show better perfomance than others and for this reason they are widely used in Finance 
[JacOll IGla03| . However, it is also known that construction of efficient Sobol’ sequences heavily 
depends on the so-called initial numbers and therefore very few Sobol’ sequence generators show 
good efficiency in practical tests [ SAKK12 ], 

Compared to Monte Carlo, Quasi Monte Carlo techniques also have some disadvantages. 
Firstly, there is no “in sample” estimation of errors: since LDS are deterministic, there is 
not a notion of probabilistic error. There have been developed some techniques, known under 
the name of randomized Quasi Monte Carlo, which introduce appropriate randomizations in 
the construction of LDS, opening up the possibility of measuring errors through a confidence 
interval while preserving the convergence rate of Quasi Monte Carlo }Gla03| . The drawback 
is the sacrifice of computational speed and, often, of some precision. Secondly, effectiveness of 
Quasi Monte Carlo depends on the integrand function, and, most importantly, the convergence 
rate can depend on the dimensionality of the problem. The latter can be seen as a big obstacle, 
since many problems in financial engineering (especially in risk management) are known to 


lr The Monte Carlo method was coined in the 1940s by John von Neumann, Stanislaw Ulam and Nicholas 
Metropolis, working on nuclear weapons (Manhattan Project) at Los Alamos National Laboratory IVN51I . 
Metropolis suggested the name Monte Carlo, referring to the Monte Carlo Casino, where Ulam’s uncle often 
gambled away his money |Met87l . Enrico Fermi is believed to have used some kind of “manual simulation” in 
the 1930s, working out numerical estimates of nuclear reactions induced by slow neutrons, with no computers 
l‘Lab66l IMet871 . 
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be high-dimensional. However, many financial applications have been reported where Quasi 
Monte Carlo outperforms standard Monte Carlo even in the presence of very high dimensions 
[PT951 [PPM ICMQ971 IKMR.Z98al IKMRZ98bl iKSHTl ISAKK12) . This fact is usually explained 
by a reduced effective dimension of the problem, with respect to its nominal dimension. The 
concept of effective dimensions was introduced in (GMQ97] , It was suggested that QMC is 
superior to MC if the effective dimension of an integrand is not too large. The notion is based 
on the ANalysis Of VAriances (ANOVA). In ILOOOl it was shown how the ANOVA components 
are linked to the effectiveness of QMC integration methods. It is important to measure the 
effective dimension in order to predict the efficiency of a Quasi Monte Carlo algorithm. Moreover, 
various techniques can be used to reduce effective dimension and, thus, improve efficiency: this 
is possible because the effective dimensiord can vary by changing the order in which the variables 
are sampled. The optimal way to achieve this can be a hard task, it could depend on the specific 
model and a general solution is not known at present. One popular choice in the financial 
literature on path-dependent option pricing [ CMQ97I. KS07] is to apply the Brownian bridge 
discretization to the simulation of the underlying stochastic process, which is based on the use 
of conditional distributions. Unlike the standard discretization, which generates values of the 
Brownian motion sequentially along the time horizon, the Brownian bridge discretization first 
generates the Brownian motion value at the terminal point, then it fills a midpoint using the value 
already found at the terminal point and then subsequent values at the successive midpoints using 
points already simulated at previous steps. In terms of QMC sampling, this simulation scheme 
means that the first coordinate of the QMC vector is used to simulate the terminal value of the 
Brownian motion, while subsequent coordinates are used to generate intermediate points. There 
are many studies which show that superior performance of the QMC approach with the Brownian 
bridge discretization in comparison with the standard discretization using MC or QMC sampling, 
in application e.g. to Asian options |jCMO~97 . [KSfl7j . However, it was pointed in |Pap01| that, 
in some cases, the Brownian bridge can perform worse than the standard discretization in QMC 
simulation. The big question is how to know with certainty which numerical scheme will provide 
superior efficiency in QMC simulation. Global Sensitivity Analysis (GSA) is the answer. 

GSA is a very powerful tool in the analysis of complex models as it offers a comprehensive 
approach to model analysis. Traditional sensitivity analysis, called local within the present 
context, applied to a function f(x ) is based on specifying a point xq in the function domain 
and then computing a derivative at x = xq. GSA instead does not require to specify a 
particular point xq in the domain, since it explores the whole domain (hence the name global). 
It also quantifies the effect of varying a given input (or set of inputs) while all other inputs are 
varied as well, providing a measure of interactions among variables. GSA is used to identify 
key parameters whose uncertainty most affects the output. This information can be used to 
rank variables, fix unessential variables and decrease problem dimensionality. Reviews of GSA 
can be found in [SKf)5aj and SAA + 10 . The variance-based method of global sensitivity indices 
developed by Sobol’ [SobOlj became very popular among practitioners due to its efficiency 
and easiness of interpretation. There are two types of Sobol’ sensitivity indices: the main 
effect indices, which estimate the individual contribution of each input parameter to the output 
variance, and the total sensitivity indices, which measure the total contribution of a single input 
factor or a group of inputs. 

For modelling and complexity reduction purposes, it is important to distinguish between 
the model nominal dimension and its effective dimension. The notions of effective dimension 
in the truncation and superposition sense were introduced in |CM097] , Further, Owen added 
the notion of “average dimension” which is more practical from the computational point of 
view I.O(K ) . Definitions and evaluations of effective dimensions are based on the knowledge of 


2 Actually, the effective dimension in the truncation sense can be reduced in this way. See Section [3] for the 
formal definition of effective dimensions. 
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Sobol’ sensitivity indices. Quite often complex mathematical models have effective dimensions 
much lower than their nominal dimensions. The knowledge of model effective dimensions is very 
important since it allows to apply various complexity reduction techniques. In the context of 
quantitative Finance, GSA can be used to estimate effective dimensions of a given problem. In 
particular, it can assess the efficiency of a particular numerical scheme (such as the Brownian 
bridge or standard discretizations). 

The paper is organized as follows: Section [2] contains a brief review on Quasi Monte Carlo 
methodology and on Low Discrepancy Sequences, with particular emphasis on financial appli¬ 
cations. Section [3] introduces GSA and the notions of effective dimensions, establishing a link 
with QMC efficiency. In Section [4] we present the results of prices and sensitivities (greeks) 
computation for selected payoffs: both GSA and convergence analysis are performed, with the 
purpose to compare MC and QMC efficiencies via a thorough error analysis. Finally, conclu¬ 
sions and directions of future work are given in Section [21 In particular we propose to apply 
our methodology to risk management issues, where a faster and smoother convergence would 
represent a great advantage in terms of both computational effort and budget. Some technical 
details are discussed in the Appendices. 

2 Monte Carlo and Quasi Monte Carlo Methods in Finance 

2.1 General motivation 

In Finance, many quantities of interest, such as prices and greeks, are defined as expectation 
values under a given probability measure, so their evaluation requires the computation of mul¬ 
tidimensional integrals of a (generally complicated) function. 

Let’s consider a generic financial instrument written on a single asset S with a single payment 
date T. We denote the instrument’s payoff at time T as V(St, 9), where St is the underlying asset 
value at time t € [0,T], and 9 is a set of relevant parameters, including instrument parameters, 
such as strikes, barriers, fixing dates of the underlying S, callable dates, payment dates, etc., 
described in the contract, and pricing parameters, such as interest rates, volatilities, correlations, 
etc., associated with the pricing model. 

Using standard no-arbitrage pricing theory, see e.g. [DufOlj . the price of the instrument at 
time t = 0 is given by 


Vo(0)=-EP[D{O,T)V(S t ,O)\F o ], 

(2.1) 

D(0, T) = exp J r(i)df^, 

(2.2) 


where (fi, J-, Q ) is a probability space with risk-neutral probability measure Q and filtration Tt 
at time t, E^[-] is the expectation with respect to Q, D(0,T) is the stochastic discount factor, 
and r(t) is the risk-neutral short spot interest rate. Notice that the values of S at intermediate 
times t before final payment date T may enter into the definition of the payoff V. 

In order to price the financial instrument, we assume a generic Wiener diffusion model for 
the dynamics of the underlying asset S, 

dS t = p(t,St)dt + a(t,S t )dW t p , (2.3) 

with initial condition So, where P is the real-world probability measure, p is the real-world 
drift, a is the volatility, and Wf is a Brownian motion under P, such that dWt ~ Zyfdt , where 
Z ~ iV(0,1) is a standard normal random variable. The solution to eq. (12.31) is given by 

S t = S 0 + [ n ( u , S u ) du+ [ a (u, S u ) dW. f , (2.4) 

Jo Jo 
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see e.g. [ Oks92] . In particular, in the Black-Scholes model the underlying asset St follows a 
simple log-normal stochastic process 


dSt = n St dt + a Sf dW t p , 


(2.5) 


with constant /i and a. The solution to equation (12.51) in a risk-neutral world (under the risk- 
neutral probability measure Q) is given b}{j 


S t = S 0 exp 



t + aW.? 


( 2 . 6 ) 


“Greeks” are derivatives of the price Vo(0) w.r.t. specihc parameters 6. They are very important 
quantities which need to be computed for hedging and risk management purposes. In the present 
work, we will consider in particular the following greeks: 


A = 

r = 


dV 0 
dS 0 ’ 
d 2 V 0 
dS$ ’ 


V = 


dVo 

da ’ 


(2.7) 

( 2 . 8 ) 

(2.9) 


called delta, gamma and vega, respectively. Notice that, in the Black-Scholes model, delta 
is exactly the hedge of the financial instrument w.r.t. the risky underlying S, and vega is a 
derivative w.r.t. a model parameter (the constant volatility a in the Black-Scholes SDE (12.51) 1. 

The solution to the pricing equation (12.11) requires the knowledge of the values of the un¬ 
derlying asset S at the relevant contract dates {Ti,... ,T n }. Such values may be obtained by 
solving the SDE (12.41) . If the SDE cannot be solved explicitly, we must resort to a discretiza¬ 
tion scheme, computing the values of S on a time grid {ti,... where t\ < t-z < ■ ■ ■ < tjj, 
and D is the number of time steps. Notice that the contract dates must be included in the 
time grid, {Ti,... ,T n } C {ti,... ,tr>}. For example, the Euler discretization scheme consists of 
approximating the integral equation (12.41) by 


Sj = Sj- 1 + fj, (tj-i,Sj-i) Atj + a (tj- 1, Sj- 1) A Wj, j = 1,... ,D 


( 2 . 10 ) 


where Atj = tj — tj- 1, A Wj = Wj — Wj- 1, and to = 0. In particular, the discretization of the 
Black-Scholes solution (|2.6D leads to 


Sj = Sj-i exp 


r ~ y ) At j + 


, 3 = 1, 


,D. 


( 2 . 11 ) 


Clearly, the price in eq. (12.11) will depend on the discretization scheme adopted. See [Kl’95 for 
the order of convergence of Euler and other discretization schemes. 

We consider two discretization schemes in eq. (12.111) : standard discretization (SD) and 
Brownian bridge discretization (BBD). In the SD scheme the Brownian motion is discretized as 
follows: 

AWj- = \fATjZj, j = 1,... ,D. (2.12) 

In the BBD scheme the first variate is used to generate the terminal value of the Brownian 
motion, while subsequent variates are used to generate intermediate points, conditioned to points 

3 We assume a constant interest rate r for simplicity. See e.g. ( BM061 . appendix B, for a generalization to 
stochastic interest rates. 
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already simulated at earlier and later time steps, according to the following formula, 


wb = o, 


Wd = \/ AtooZ] , 


Wi = Wi + W k + 


At k jAt 




At 


ki 


At 


ki 




Zu 


ti < tj < t k , l — 2,... ,D, 


(2.13) 


where At a & = t a — t),. Unlike the SD scheme, which generates the Brownian motion sequentially 
across time steps, the BBD scheme uses different orderings: as a result, the variance in the 
stochastic part of (12.1311 is smaller than that in (12.121) for the same time steps, so that the first 
few points contain most of the variance. Both schemes have the same variance, hence their MC 
convergence rates are the same, but QMC sampling shows different efficiencies for SD and BBD, 
which will be discussed in the following sections. 

The number D of time steps required in the discretization of the SDE (12.1011 is the nominal 
dimension of the computational problem: indeed, the expectation value in (12.111 is formally an 
integral of the payoff, regarded as a function of Z\,... , Zy). In general, financial instruments may 
depend on multiple underlying assets S' 1 ,... ,S AI : in this case, the dimension of the problem 
is given by D x M. In conclusion, the pricing problem (12.111 is reduced to the evaluation of 
high-dimensional integrals. This motivates the use of Monte Carlo techniques. 

Throughout this work, we will focus on the relative effects of the dimension D and of the 
discretization schemes on the MC and QMC simulations. Thus, we will assume a simple Black- 
Scholes underlying dynamics for simplicity. This choice will be also useful as a reference case 
to interpret further results based on more complex dynamic^. We stress that using simple and 
solvable dynamics is an approximation often used in risk management practice for risk measures 
calculation on large portfolios with multiple underlying risk factors, because of computational 
bottlenecks. 


2.2 Pseudo Random Numbers and Low Discrepancy Sequences 


Standard gaussian numbers Zj are computed using transformation of uniform variates xj ~ 
i.i.d. U( 0,1), 

Zj = j = l,...,D, (2.14) 

where T -1 is the inverse cumulative distribution function of the standard normal distribution. 
Hence the pricing problem (12.111 can be reduced to the evaluation of integrals of the following 
generic form 

V= [ f(x)d D x, (2.15) 

Jh d 

where H D = [0, 1} D is the D —dimensional unit hypercube. The standard Monte Carlo estimator 
of (12.1511 has the form 

1 N 

Vn - (2.16) 
V k =1 

where {x k }^ =1 is a sequence of N random points in H D . Sequences {x k } k=1 are produced 
by appropriate Random Number Generators (RNGs). In particular, Pseudo Random Number 
Generators (PRNGs) are computer algorithms that produce deterministic sequences of pseudo 
random numbers (PRNs) mimicking the properties of true random sequences. Such sequences are 
completely determined by a set of initial values, called the PRNG’s state. Thus, pseudo random 
sequences are reproducible, using the same set of state variables. PRNGs are characterized 

4 For example, we could introduce jumps or Heston dynamics, see e.g. IWilOBI . 
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by the seed, i.e. a random number used to initialize the PRNG, the period, i.e. the maximum 
length, over all possible state variables, of the sequence without repetition, and the distribution 
of the generated random numbers, which is generally uniform [0,1). The most famous PRNG is 
the Mersenne Twister [MN98j . with the longest period of 2 19937 — 1 and good equidistribution 
properties guaranteed up to, at least, 623 dimensions. Pseudo random sequences are known 
to be plagued by clustering: since new points are added randomly, they don’t necessarily fill 
the gaps among previously sampled points. This fact causes a rather slow convergence rate. 
Consider an integration error 

s=\V - V N \. (2.17) 


By the Central Limit Theorem the root mean square error of the Monte Carlo method is 

_2mV 2 _ a f 


s MC = = 


VN' 


(2.18) 


where <7/ is the standard deviation of f(x). Although emc does not depend on the dimension 
D, as in the case of lattice integration on a regular grid, it decreases slowly with increasing 
N. Variance reduction techniques, such as antithetic variables [JacOll IGla03| , only affect the 
numerator in (12.181) . 

In order to increase the rate of convergence, that is to increase the power of N in the 
denominator of (12.181) . one has to resort to Low Discrepancy Sequences (LDS), also called Quasi 
Random Numbers (QRNs), instead of PRNs. The discrepancy of a sequence {xk}^ = i is a 
measure of how inhomogeneously the sequence is distributed inside the unit hypercube H D . 
Formally, it is defined by [Jac nu 


T>n{x\, • • -,x N ) = sup 


n 


[S d (£),x 1 ,...,x n ] 


N 


~ m(0 


D 


(2.19) 


S D (0 = [0,6) x ■ ■ ■ x [0,6,) C H d , m(() = I] 6, 

j =i 


where 


n 


N N D 

[*S D (0, *i, ■ ■ ■, *jv] = = S II 

k =1 k=l j=l 


( 2 . 20 ) 


is the number of sampled points that are contained in hyper-rectangle S D C H D . It can be shown 
that the expected discrepancy of a pseudo random sequence is of the order of ln(ln N)/\/N. A 
Low Discrepancy Sequence is a sequence {x^}^ =1 in H D such that, for any N > 1, the first N 
points Xi ,..., xn satisfy inequality 

V%{xi,...,x N ) < c(D) ln / , (2.21) 

for some constant c(D) depending only on D |Nie88j . Unlike PRNGs, Low Discrepancy Se¬ 
quences are deterministic sets of points. They are typically constructed using number theoretical 
methods. They are designed to cover the unit hypercube as uniformly as possible. In the case 
of sequential sampling, new points take into account the positions of already sampled points 
and fill the gaps between them. Notice that a regular grid of points in H D does not ensure low 
discrepancy, since projecting adjacent dimensions easily produces overlapping points. 

A Quasi Monte Carlo (QMC) estimator of the integral (12.151) is of the form (12.161) with the 
only difference that the sequence is sampled using LDS instead of PRNs. An upper 

bound for the QMC integration error is given by the Koksma-Hlawka inequality 

SQMC < V(f)V% = O (^) 


( 2 . 22 ) 

























where V (/) is the variation of the integrand function in the sense of Hardy and Krause, which 
is finite for functions of bounded variation [KFSMl lJ. The convergence rate of ([2.22D is asymp¬ 
totically faster than 112.181) . but it is rather slow for feasible N. Moreover, it depends on the 
dimensionality D. However, eq. (|2.22l) is just an upper bound: what is observed in most 
numerical tests [ KFSM1L CMQ97] is a power law 

eqmc ~ , (2.23) 

where the value of a depends on the model function and, therefore, is not a priori determined 
as for MC. When a > 0.5 the QMC method outperforms standard MC: this situation turns out 
to be quite common in financial problems. We will measure a for some representative financial 
instruments, showing that its value can be very close to 1 when the effective dimension of / 
is low, irrespective of the nominal dimension D. The concept of effective dimension, and the 
methodology to compute it, will be introduced in the following sections. 

We stress that, since LDS are deterministic, there are no statistical measures like variances 
associated with them. Hence, the constant c in l[2.23l) is not a variance and l[2.23l) does not have 
a probabilistic interpretation as for standard MC. To overcome this limitation, Owen suggested 
to introduce randomization into LDS at the same time preserving their superiority to PRN 
uniformity properties |Owe93j . Such LDS became known as scrambled (see also |Gla03| ). In 
practice, the integration error for both MC and QMC methods for any fixed N can be estimated 
by computing the following error averaged over L independent runs: 


£N = 



E(r-0 2 . 


(2.24) 


where V is the exact, or estimated at a very large extreme value of IV —>• oo, value of the integral 
and Vfi is the simulated value for the Hh run, performed using N PRNs, LDS, or scrambled 
LDS. For MC and QMC based on scrambled LDS, runs based on different seed points are 
statistically independent. In the case of QMC, different runs are obtained using non overlapping 
sections of the LDS. Actually, scrambling LDSs weakens the smoothness and stability properties 
of the Monte Carlo convergence, as we will see in Section 14.51 Hence, in this paper we will use 
the approach based on non-overlapping LDSs, as in [SK05b| . 

The most known LDS are Halton, Faure, Niederreiter and Sobol’ sequences. Sobol’ se¬ 
quences, also called LPt sequences or (t, s) sequences in base 2 [Nie88 j. became the most known 
and widely used LDS in finance due to their efficiency [JacOll !Gla03j . Sobol’ sequences were 
constructed under the following requirements fSob67j : 

1. Best uniformity of distribution as N —>• oo. 


2. Good distribution for fairly small initial sets. 

3. A very fast computational algorithm. 

The efficiency of Sobol’ LDS depend on the so-called initialisation numbers. In this work we 
used SobolSeq8192 generator provided by BRODA [BROj . SobolSeq is an implementation of 
the 8192 dimensional Sobol’ sequences with modified initialisation numbers. Sobol’ sequences 
produced by SobolSeq8192 can be up to and including dimension 2 13 , and satisfy additional uni¬ 
formity properties: Property A for all dimensions and Property A’ for adjacent dimensions (see 
SAKK 12 for detail^). It has been found in [SAKK 12) that SobolSeq generator outperforms 
all other known LDS generators both in speed and accuracy. 

5 BRODA releases also SobolSeq32000 and SobolSeq64000 
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3 Global Sensitivity Analysis and Effective Dimensions 


As we mentioned in the Introduction and Section m, effective dimension is the key to explain 
the superior efficiency of QMC w.r.t. MC. Hence, it is crucial to develop techniques to estimate 
the effective dimension and to find the most important variables in a MC simulation. 

The variance-based method of global sensitivity indices developed by Sobol’ became very 
popular among practitioners due to its efficiency and easiness of interpretation [SK05al ISAA + 10 . 
There are two types of Sobol’ sensitivity indices: the main effect indices, which estimate the 
individual contribution of each input parameter to the output variance, and the total sensitivity 
indices, which measure the total contribution of a single input factor or a group of inputs. Sobol’ 
indices can be used to rank variables in order of importance, to identify non-important variables, 
which can then be fixed at their nominal values to reduce model complexity, and to analyze the 
efficiency of various numerical schemes. 

Consider a mathematical model described by an integrable function f(x), where the input 
x = (aji,..., xd) is taken in a D-dimensional domain 17 and the output is a scalar. Without loss 
of generality, we choose 17 to be the unit hypercube H D . The input variables x\,... ,xd can, 
then, be regarded as independent uniform random variables each defined in the unit interval 
[0,1]. The starting point of global sensitivity analysis (GSA) is the analysis of variance (ANOVA) 
decomposition of the model function, 


f{x) = fo + ^2fi{xi) + + • • • + / i 2 - d ( zi , ...,x d ). ( 3 . 1 ) 

i i<j 


The expansion (13.11) is unique, provided that 


f 


fil-.-is { x ii J • • ■ ; x i a )dxi k =0, Mk= l,...,s. 


(3.2) 


The ANOVA decomposition expands the function / into a sum of terms, each depending on 
an increasing number of variables: a generic component fi 1 —i a (xi 1 ,...,Xi a ), depending on s 
variables, is called an s-order term. It follows from (13.21) that the ANOVA decomposition is 
orthogonal and that its terms can be explicitly found as follows, 


fo = [ f{x)d D x, 
Jh d 

•(*<)=/ /wn 

JHv - 1 


dx k /o> 


(3.3) 


fij(xi,Xj ) = / f(x) TT dx k - fo - fi(xi ) - fj(xj), 

Jhd ~ 2 kA 


and so on. If / is square-integrable, its variance decomposes into a sum of partial variances: 


^2 = 


where 


U; 


^ + + • • • + 
i i<j 

[ fi 1 ...i 3 (x il ,...,x is )d. 

■Jo 


cr 12 ■■■£>> 


l i ■ ■ ■ dxi s . 


Sobol’ sensitivity indices are defined as 


C _ T 1 

°h-i s - 2 


(T- 


(3.4) 


(3.5) 


(3.6) 
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and measure the fraction of total variance accounted by each fi 1 ...i s term of the ANOVA decom¬ 
position. From (13.41) it follows that all Sobol’ indices are non negative and normalized to 1. First 
order Sobol’ indices S t measure the effect of single variables Xi on the output function; second 
order Sobol’ indices Sij measure the interactions between pairs of variables, i.e. the fraction 
of total variance due to variables Xi and Xj which cannot be explained by a sum of effects of 
single variables; higher order Sobol’ indices Sq...^, with s > 2, measure the interactions among 
multiple variables, i.e. the fraction of total variance due to variables x^ ,... , Xj s which cannot 
be explained by a sum of effects of single variables or lower order interactions. 

The calculation of Sobol’ sensitivity indices in eq. (13.61) requires, in principle, 2 D valuations 
of the multi-dimensional integrals in eq. (|3.5D . which is a very cumbersome, or even impossible, 
computational task. Furthermore, for practical purposes, and in particular when the function 
/ has low order interactions, it is not actually necessary to know all the possible Sobol’ indices, 
but just an appropriate selection of them. Thus, it is very useful to introduce Sobol’ indices for 
subsets of variables and total Sobol’ indices. Let y = {xq,..., Xj m } C x, 1 < i\ < ..., < i m < D, 
be a subset of x, and z = y^ C x its complementary subset, and define 


D 


s » = £ £ 

s=l (h<—<i s )eK 




QtOt -| Q 

Oy - 1 D Z , 


(3.7) 


where K = {i\ ,..., i m }. Notice that 0 < S y < Sy 0t < 1. The quantity Sy 0t — S y accounts for all 
the interactions between the variables in subsets y and z. It turns out that there exist efficient 
formulas which allow to avoid the knowledge of ANOVA components and to compute Sobol’ 
indices directly from the values of function / [SobOl] . These formulas are based on the following 
integrals, 


Sy ~G 2 


[. f(y', A - fo\ [f(yz) - f{y , z)\dy dz dy'dz 1 , 

Sy>t J 0 ~ Z ^ dy dz dy ' ’ 

A = [ f 2 {y, z)dy dz — /o , 

Jo 

fo= f(y,z)dydz, 

Jo 


(3.8) 


where the integration variables are the components of the vectors y, z, y', z', such that x = yL) z, 
and the first two integrals depend on the choice of y. Such integrals can be evaluated, in general, 
via MC/QMC techniques [KFSMlll iSal02] . 

Furthermore, usually enough information is already given by the first order indices Si and 
by corresponding total effect indices Sj ot , linked to a single variable y = {x^}. For these Sobol’ 
indices, it’s easy to see that 

• Sj ot = 0: the output function does not depend on x*, 


• Si = 1: the output function depends only on Xj, 

• Si = Sf ot : there is no interaction between Xj and other variables. 


Notice that just D + 2 function evaluations for each MC trial are necessary to compute all S) and 
S\ ot indices in eqs. (13.81) : one function evaluation at point x = {y, z}, one at point x' = {y' , z'}, 
and D evaluations at points x" = {y\ z} , V y 1 = {xj} , i = 1,..., D. 
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We stress that the approach presented above is applicable only to the case of independent 
input variables, which admits a unique ANOVA decomposition. In the case of dependent (cor¬ 
related) input variables, the computation of variance-based global sensitivity indices is more 
involved. A generalization of GSA to dependent variables can be found in [KTA12j . 

We finally come to the notion of effective dimensions, firstly introduced in [CM097] , Let \y\ 
be the cardinality of a subset of variables y. The effective dimension in the superposition sense, 
for a function / of D variables, is the smallest integer ds such that 

E S y >l-e ( 3 - 9 ) 

0<|y|<d s 

for some threshold e (arbitrary and usually chosen to be less than 5%). If a function has an 
effective dimension dg in the superposition sense, it can be approximated by a sum of ds~ 
dimensional terms, with an approximation error below e. 

The effective dimension in the truncation sense is the smallest integer dx such that 

E S y > 1-e. (3.10) 

yC{l,2,...,dr} 

The effective dimension ds does not depend on the order of sampling of variables, while dx does. 
In general, the following inequality holds, 


ds S dx £ D. 


(3.11) 


Effective dimensions can be estimated solely from indices S) and Sj ot using eqs. (13.81) with 
y = i, as described in [KFSMlTj , where relationships among such indices are used to classify 
functions in three categories according to their dependence on variables. For the so-called type 
A functions, variables are not all equally important and the effective dimension in the truncation 
sense dx is small, such that dg < dx D. They are characterized by the following relationship 


S.t ot 5 

v > - 


'tot 


\y\ 


(3.12) 


where y C x is a leading subset of variables, z = y^ C x its complementary subset. Functions 
with equally important variables have dx — D and they can be further divided in two groups: 
type B and C functions. Type B functions have dominant low-order interactions and small 
effective dimension in the superposition sense ds, such that ds <C dx — D. For such functions, 
Sobol’ indices satisfy the following relationships: 


D 

Si ~ Sj ot , Vi = l,...,D, E 5 *- 1 ' ( 3 ‘ 13 ) 

2—1 

Type C functions have dominant higher-order interactions 

D 

S t < Sf 4 , E < :l ( 3 - 14 ) 

2 — 1 

and effective dimensions ds — dx — D. This classification is summarized in Table [T| 

Owen introduced in |Qwe03| the notion of the average dimension d^, which can assume 
fractional values, defined as 

d A -= E \y\ s v’ ( 3 - 15 ) 

0<\y\<D 
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Type 

Description 

Relationship between SI 

Eff. dimensions 

A 

Few important variables 

ST/\y \» sr/M 

ds < dr "C D 

B 

Low-order interactions 


ds dr — D 

C 

High-order interactions 

S z ~ Sj,Si <C S\ ot ,Vi,j 

ds — dx — D 


Table 1: Classification of functions w.r.t. their dependence on variables, based on GSA. 


and showed that it can be rather straightforwardly computed as 

D 

d A = Y J S? t . (3.16) 

%— 1 

It has been suggested in |SS14] that QMC should outperform MC when d A < 3. This is confirmed 
in our findings, see Section fOl 

It has been proved in many works [KFSM11 , C.MQ97. Owef)3j that QMC outperforms MC 
regardless of the nominal dimension whenever the effective dimension is low in one or more 
senses. Hence, in the case of type A and type B functions (we assume that functions are 
sufficiently smooth), QMC always outperform MC, while for type C functions the two methods 
are expected to have similar efficiency. Actually, type A and B functions are very common 
in financial problems. We also note that the performance of the QMC method for Type A 
functions sometimes, but not always, can be greatly improved by using effective dimension 
reduction techniques, such as Brownian bridge, which will be demonstrated in the following 
section. 


4 Test Cases and Numerical Results 

In this section we apply MC and QMC techniques to high-dimensional pricing problems. Our 
aim is to test the efficiency of QMC with respect to standard MC in computing prices and 
greeks (delta, gamma, vega) for selected payoffs V with increasing degree of complexity and 
path-dependency. 


4.1 Selected Payoffs and Test Set-Up 

We selected the following instruments as test cases. 

1. European call: 


V = max(S , £) — K, 0). 


2. Asian call: 


1 ID 


V = max(5 — K, 0), S = ‘-b 

3. Double knock-out: 

V = rna x(S D - K, 0) l{R l <s. j <B u }, V? = 1,..., D. 

4. Cliquet: 


D 


V = max < max 
3 =i 


0, min ( C, 


Sj ~ Si-i 


Si-1 


,F 


(4.1) 


(4.2) 


(4.3) 


(4.4) 
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In the above definitions, K denotes the strike price, Bi and B u are the values of the lower and 
upper barrier, respectively, C is a local cap and F is a global floor. In all test cases we use the 
following payoff parameters: 

• maturity: T = 1, 

• strike: K = 100, 

• lower barrier: Bi = 0.5 So, 

• upper barrier: B u = 1.5 So, 

• global floor: F = 0.16, 

• local cap: C = 0.08. 

Such selection guarantees an increasing level of complexity and path-dependency. The European 
call is included just as a simple reference case, for which analytical formulas are available for 
price and greeks, see e.g. |Wil06| . The Asian call with arithmetic average is the simplest and 
most diffused non-European payoff; we choose geometric average payoff such that analytical 
formulas are availably. The double barrier is another very diffused payoff with stronger path- 
dependency. Finally, the Cliquet option is a typical strongly path-dependent payoff based on 
the performance of the underlying stock. Clearly, many other possible payoffs could be added to 
the test (e.g. autocallable), but we think that such selection should be complete enough to cover 
most of the path-dependency characteristics relevant in the Monte Carlo simulation. We assume 
that the underlying process St follows a geometric Brownian motion as described in Section 12.11 
with the following model parameters: 

• spot: So = 100, 

• volatility: a = 0.3, 

• number of time steps: D = 32. 

The process St is discretized across D time steps {t\ < ■ ■ ■ < tj < ■ ■ ■ < to}, so that So is its 
value at maturity. Recall that, in the single asset case, the number of time simulation steps is 
equal to the dimension of the path-dependent simulation. As discussed at the end of Section O 
we choose a simple dynamics for St because our main goal to compare MC and QMC simulations 
w.r.t. the effect of the dimension D and of the discretization schemes. 

The numerical computations are performed in Mat lab using three different sampling tech¬ 
niques: 

• MC+SD + antithetic variables + Mersenne Twister generator, 

• QMC+SD + SobolSeq8192 generator, 

• QMC+BBD + SobolSeq8192 generator. 

The notations for the simulation parameters are: 

• N: number of simulated paths for the underlying, 

• D: number of time steps used to discretize each underlying’s path, 

• L : number of independent runs. 

6 See e.g. IWil06l and references therein. 
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Notice that, using the Black-Scholes model, the number D of time steps is also the nominal 
dimension of the MC simulation. Following the specifics of Sobol’ sequences, we take N = 2 P , 
where p is an integer, since this guarantees the lowest discrepancy properties. 

Simulation errors £jv are analyzed by computing the root mean square error (RMSE) as 
defined by (12.24k where V is a reference value of prices or greeks given by analytical formulas 
(for European and geometric Asian options) or simulated with a large number of scenarios 
(N = 2 23 ) (for Double Knock-out and Cliquet options). To assess and compare performance 
of MC and QMC methods with different discretization schemes, we compute the scaling of the 
RMSE as a function of N by fitting the function with a power law cN~ a (12.2311 . In the MC 
case, the value of a is expected to be 0.5 in all situations, while in the QMC case it is expected 
to be higher than 0.5 for Type A and B functions. 

Finally, greeks for the payoffs above are computed via finite differences, using central differ¬ 
ence formulas for delta, gamma and vega, with shift parameter e, 


dvf V 0 v (Sq + h)~ Vf(S 0 -h) 

F dS 0 - 2 h 

d 2 Vf Vf (S 0 + h) - 2Vf (So) + V 0 v (So - h ) 

v ds 2 ~ v 

v d_Yl^ V£(o + h)-V£(a-h) 
da 2 h ’ 

where the increment h is chosen to be h = eSo, f° r delta and gamma, and h = e, for vega, for 
a given “shift parameter” e. Notice that the calculation of price and three greeks using eqs. 
(USD and payoffs (14.1114.41) above requires N p = 5 + 5 + 5 + 3 = 18 functions evaluations (the 
Cliquet has null delta and gamma). In the MC simulations for greeks we use path recycling of 
both pseudo random and LDS sequences to minimize the variance of the greeks, as suggested 
in [Jac m and m a03j . Notice that the analysis of the RMSE for greeks is, in general, more 
complex than that for prices, since the variance of the MC simulation mixes with the bias due 
to the approximation of derivatives with finite differences with shift e. We discuss how to deal 
with this issue in Appendix [XI 


4.2 Global Sensitivity Analysis for Prices and Greeks 

Sobol’ indices S t and Sj ot are computed for both the standard and Brownian bridge dis¬ 
cretizations using eqs. (13.81) . where / is the relevant model function (the instrument payoff 
or a greek with finite differences) and y = {xj}, y' = {x'}, z = {xi,..., Xj_i, ajj+i,..., xd}, 
z' = {x^,..., x\_ x , x' +1 ,..., x'h]. Here x* are the uniform variates x* ~ i.i.d. U [0,1] used in 
eq. (12.141) . The integrals in eqs. (13.81) are computed using QMC simulation with the following 
parameters: 

• number of simulations: N = 2 , 

• shift parameter for finite difference^ e = 10 —4 ,10 -3 ,10~ 2 . 

Effective dimensions are estimated in the following way: 

• The effective dimension in the truncation sense dx is computed using inequality (13.121) . 
looking for a minimal set of variables y = {xi,..., Xd T } such that the quantity S^\y\/Sy 0t |^| 
is smaller than 1%. Since the calculation of dx depends on the order of sampling variables, 
the result depends on the discretization scheme used, that is SD or BBD. 

7 See Appendix E 
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Payoff 

Function 

Si/Sl ot 

Si Si 

dx 

ds 

dA 

Effect of e 

European 

Price 

0.49 

0.68 

32 

< 32 

1.40 

- 


Delta 

0.26—>0.23 

0.77 

32 

< 32 

3.2 

small 


Gamma 

10" 4 -> 10~ 2 

10" 4 -> 10" 2 

32 

32 

32 

high 


Vega 

0.33 

0.543 

32 

< 32 

1.64 

no 

Asian 

Price 

0.54—>0.43 

0.714 

< 32 

< 32 

1.38 

- 


Delta 

0.32—> 10 -2 

0.71—>0.74 

32 

< 32 

3.5 

small 


Gamma 

10" 4 -> 10 -2 

10 -4 -> 10" 2 

32 

32 

31 -> 25 

high 


Vega 

0.42—>0.01 

0.611 

< 32 

< 32 

1.57 

no 

Double KO 

Price 

0.01—>0.15 

0.22 

32 

< 32 

8.5 

- 


Delta 

0.01-> 0.12 

0.22 

32 

< 32 

7.6 

no 


Gamma 

10" 5 -> 10“ 7 

10" 4 -> 10" 2 

32 

32 

31.2 -> 29.8 

high 


Vega 

10" 5 -> 10“ 8 

10~ 4 -> 10" 2 

32 

32 

28 

high 

Cliquet 

Price 

1 

1 

32 

1 

1 

- 


Vega 

1 

1 

32 

1 

1 

no 


Table 2: Summary of GSA metrics and effective dimensions of prices and greeks for SD scheme. Arrow 
“—>” in the column for Si/S '| ot denotes the change in the value with the increase of index i and/or 
with the increase of shift parameter e; in the column for Y/i <5* it denotes the change in the value 
with the increase of shift parameter e. The numerical computation of the figures in this table required 
N p x (D + 2) x N e x TV = 18 x 34 x 3 x 2 17 = 240,648,192 function evaluations. We show significant 
digits only, we do not show MC errors because of limited space. 


• The effective dimension in the superposition sense ds is estimated using dimension dx as 
an upper bound according to inequality (13.111) . In order to distinguish between Type B 
and Type C functions, we look at ratios Si/S\ ot and Y2i according to eqs. (13.131) . (13.141) . 

• The effective average dimension dA is computed according to eq. (13.161) . 

The results of GSA for the SD are shown in Figures [TB Measures based on Sobol’ indices are 
provided in Table [2j These measures are used to compute effective dimensions and to classify 
integrands in (12.151) corresponding to prices and greeks according to Table [0 
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(a) Price 


(b) Delta 
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0.03 
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0.015 
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35 


(c) Gamma (d) Vega 

Figure 1: European call option price (a) and greeks (6), (c), (d), SD, D = 32. First order Sobol’ indices 
Si and total sensitivity indices Sl ot versus time step i. 
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(c) Gamma 

Figure 2: Asian call option. 



(b) Delta 
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(a) Price 


(b) Delta 




(c) Gamma 


(d) Vega 


Figure 3: 


Double Knock-out call option. Parameters as in Figure [Tj 



(a) Price 


(b) Vega 


Figure 4: Cliquet option. 


Parameters as in Figure [lj Delta and gamma are null for Cliquet options. 
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From these results we draw the following conclusions. 

1. European option (Figure [I]): price, delta and vega are type B functions, while gamma is 
type C function. 

2. Asian option (Figure [2]): price, and vega are type B functions, while delta and gamma are 
type C function. 

3. Double KO option (Figure [3]) : price and all greeks are type C functions. 

4. Cliquet option (Figure 0]): price and vega are type B functions with ds = 1 (delta and 
gamma for a Cliquet option are null). We recall that ds = 1 means that there are no 
interactions among variables. 

The analogous results of GSA for BBD are shown in Figures [5]|8] and in Tabled 


Payoff 

Function 

Si/Sl ot 

EiSi 

dx 

ds 

dA 

Effect of e 

European 

Price 

1 

i 

i 

1 

1 

- 


Delta 

1 

i 

i 

1 

1 

no 


Gamma 

1 

i 

i 

1 

1 

no 


Vega 

1 

i 

i 

1 

1 

no 

Asian 

Price 

0.853—>0.4 

0.875 

2 

<2 

1.13 

- 


Delta 

0.733—>0.01 

0.778 

4 

<4 

1.68 -> 1.43 

small 


Gamma 

10" 2 -> 10 -4 

0.022 —> 10" 4 

32 

32 

31 -> 8 

high 


Vega 

0.802—>0.03 

0.827 

2 

<2 

1.20 

no 

Double KO 

Price 

0.70—>0.01 

0.70 

~ 2 

<2 

1.63 

- 


Delta 

0.83—>0.01 

0.83 

2 

<2 

1.37 

no 


Gamma 

1 

1 -> 0.95 

1 

1 

1.0 

small 


Vega 

10~ 4 -> 0.2 

10" 6 -> 10" 4 

32 

32 

4.8 -> 3.9 

high 

Cliquet 

Price 

0.978—>0.2 

0.892 

~ 2 

< 2 

1.19 

- 


Vega 

0.595—>0.001 

0.32 

~ 32 

< 32 

2.6 

no 


Table 3: Summary of GSA metrics and effective dimensions of prices and greeks for BBD scheme. 
Details as in Table [2] 


20 








Figure 5 

indices Si 



(c) Gamma (d) Vega 

: European call option price (a) and greeks (b),(c),(d), BBD, D = 32. First order Sobol’ 
and total sensitivity indices Sl ot versus time step i. 


21 



























(b) Delta 



(c) Gamma (d) Vega 

Figure 6: Asian call option. Details as in Figure [5] 
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(a) Price 

1.2 *- 1 - 1 - 1 - 1 - 1 - 1 - 

s with e = 10~ 4 
:ot with e = 10“ 4 
withe= 1CT 3 
:ot with e = 10“ 3 
with£= 10“ 2 
01 with £= 10“ 2 



-0.2- 1 - 1 - 1 - 1 - 1 - 1 - 

0 5 10 15 20 25 30 35 



(b) Delta 




(c) Gamma 


(d) Vega 


Figure 7: Double Knock-out call option. Details as in Figure [5] 




(a) Price (b) Vega 

Figure 8: Cliquet option. Details as in Figure [5] 


23 











































From these results we draw the following conclusions. 

1. European option (Figure [5]): price and all greeks are type A functions with ds = 1. The 
value of sensitivity indexes for the first input, corresponding to the terminal value t = T, 
is ~ 1, while the following variables have sensitivity indexes ~ 0. Clearly, BBD is much 
more efficient than SD. 

2. Asian option (Figure [6]): price, delta and vega are type A functions. Comments as for the 
European option above. Gamma remains a type C function as for SD. 

3. Double KO option (Figure 0: price, delta and gamma are type A functions. Comments 
as for the Asian option above. Vega remains a type C function as for SD. 

4. Cliquet option (Figure 0: price is a type A function. Similarly to the European option, 
the value of sensitivity indexes for the first input, corresponding to the terminal value 
t = T, is ~ 1, while the following values of Si are ~ 0. Vega is a type C function, since the 
ratio Si/S\ ot reaches small values revealing interacting variables. Thus in this case BBD 
is much less efficient than SD. 

In conclusion, prices and greeks are always Type B or C functions for QMC+SD (Table [2]), while 
they are predominantly Type A functions, with a few exceptions, for QMC+BBD (Table [3]). In 
most cases switching from SD to BBD reduces the effective dimension in the truncation sense 
dT- 

The different efficiency of QMC+BBD vs QMC+SD is completely explained by the properties 
of Sobol’ low discrepancy sequences. The initial coordinates of Sobol’ LDS are much better 
distributed than the later high dimensional coordinates [Gla03[ [CM097]. The BBD changes 
the order in which inputs (linked with time steps) are sampled. As follows from GSA, in most 
cases for BBD the low index variables (terminal values of time steps, mid-values and so on) are 
much more important than higher index variables. The BBD uses lower index, well distributed 
coordinates from each D-dimensional LDS vector to determine most of the structure of a path, 
and reserves other coordinates to fill in finer details. That is, well distributed coordinates are 
used for important variables and other not so well distributed coordinates are used for far less 
important variables. This results in a significantly improved accuracy of QMC integration. 
However, this technique does not always improve the efficiency of the QMC method as e.g. for 
Cliquet options: in this case GSA reveals that for SD all inputs are equally important and, 
moreover, there are no interactions among them, which is an ideal case for application of Sobol’ 
low discrepancy sequences; the BBD, on the other hand, favoring higher index variables destroys 
independence of inputs introducing interactions, which leads to higher values of ds and dA- As 
a result, we observe degradation in performance of the QMC method. 

4.3 Performance Analysis 

In this section we compare the relative performances of MC and QMC techniques. This analysis 
is crucial to establish if QMC outperforms MC, and in what sense. 

Firstly, following the suggestion of [Jac m, Section 14.4, we analyze convergence diagrams 
for prices and greeks, showing the dependence of the MC simulation error upon the number of 
MC paths. The results for the four payoffs are shown in Figures 1911121. 
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(a) Price (b) Delta 




(c) Gamma 


(d) Vega 


Figure 9: European call option price (a) and greeks (b), (c), (d) convergence diagrams versus number of 
simulated paths for MC+SD with antithetic variables (solid blue line), QMC+SD (solid green line) and 
QMC+BBD (solid red line). Shaded areas represent 3-sigma errors around the corresponding run (solid 
line). 1% and 0.1% accuracy regions are marked by horizontal black solid and dashed lines, respectively. 
Number of dimensions is D = 32. Shift parameter is e = 10~ 3 . 
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(a) Price 



(c) Gamma 



(b) Delta 



(d) Vega 


Figure 10: Asian call option, e = 5 x 10 3 . Other details as in Figure[9l 
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(a) Price (b) Delta 




(c) Gamma (d) Vega 

Figure 11: Double Knock-out call option. Details as in Figure [TO] 




(a) Price (b) Vega 

Figure 12: Cliquet option. Details as in Figure [lU] 
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We observe what follows. 


1. European option (Figure [9]): QMC+BBD outperforms both QMC+SD and MC+SD in all 
cases (the 3-sigma regions for QMC+BBD are systematically smaller). We also note that 
for price, delta and vega, the QMC+BBD convergence is practically monotonic, which 
makes on-line error approximation possible. For gamma, the QMC+BBD convergence is 
much less oscillating than for MC+SD. 

2. Asian option (Figure [TUI): QMC+BBD outperforms both QMC+SD and MC+SD for price 
and vega. For delta, QMC with both SD and BBD is marginally better than MC+SD. 
For gamma, QMC with both SD and BBD has nearly the same efficiency as MC+SD. The 
QMC+BBD convergence is also smoother for price, delta and vega. 

3. Double KO option (Figure [TTj) : QMC+BBD outperforms both QMC+SD and MC+SD in 
all cases. 

4. Cliquet option (Figure fT2l): QMC+SD outperforms QMC+BBD and MC+SD in all cases. 
QMC+BBD outperforms MC+SD only for price. 

Next, we analyze the relative performance of QMC vs MC in terms of convergence rate. 
We plot in Figures fT3lfl6l the root mean square error, eq. (12.241) . versus the number of MC 
scenarios N in Log-Log scale. In all our tests we have chosen an appropriate range for N such 
that, in the computation of greeks, the bias term is negligible with respect to the variance term 
(see Appendix [A] for details). Hence, the observed relations are, with good accuracy, linear, 
therefore the power law (12.231) is confirmed, and the convergence rates a can be extracted as the 
slopes of the regression lines. Furthermore, also the intercepts of regression lines provide useful 
information about the efficiency of the QMC and MC methods: in fact, lower intercepts mean 
that the simulated value starts closer to the exact value. The resulting slopes and intercepts 
from linear regression are presented in tabs. 0] and [5] for all test cases. 


Payoff 

Function 

MC+SD 

QMC+SD 

QMC+BBD 

European 

Price 

- 0.3 ± 0.1 

- 0.2 ± 0.1 

- 0.84 ± 0.01 


Delta 

- 2.0 ± 0.1 

- 1.6 ± 0.1 

- 2.64 ± 0.01 


Gamma 

- 2.0 ± 0.1 

- 2.1 ± 0.1 

- 2.8 ± 0.2 


Vega 

0.4 ± 0.1 

0.4 ± 0.1 

- 0.13 ± 0.01 

Asian 

Price 

- 0.5 ± 0.1 

- 0.5 ± 0.1 

- 1.0 ± 0.1 


Delta 

- 2.2 ± 0.1 

- 1.7 ± 0.1 

- 1.9 ± 0.1 


Gamma 

- 2.1 ± 0.2 

- 2.1 ± 0.1 

- 2.0 ± 0.1 


Vega 

0.1 ± 0.1 

- 0.1 ± 0.1 

- 0.4 ± 0.1 

Double KO 

Price 

- 0.4 ± 0.1 

- 0.3 ± 0.1 

- 0.7 ± 0.1 


Delta 

- 1.8 ± 0.1 

- 1.6 ± 0.1 

- 2.1 ± 0.1 


Gamma 

- 2.4 ± 0.1 

2.1 ± 0.1 

- 2.9 ± 0.1 


Vega 

1.1 ± 0.1 

1.3 ± 0.1 

1.3 ± 0.2 

Cliquet 

Price 

- 2.4 ± 0.1 

- 3.2 ± 0.1 

- 2.5 ± 0.3 


Vega 

- 2.0 ± 0.1 

- 2.7 ± 0.1 

- 1.7 ± 0.2 


Table 4: Intercepts from linear regression with their errors, for MC+SD with antithetic variables, 
QMC+SD and QMC+BBD, L = 30 runs. Results are shown for N = 10 2,5 paths. 
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Payoff 

Function 

MC+SD 

QMC+SD 

QMC+BBD 

European 

Price 

Delta 

Gamma 

Vega 

-0.46 ± 0.03 
-0.49 ± 0.03 
-0.51 ± 0.02 
-0.45 ± 0.03 

-0.71 ± 0.03 
-0.56 ± 0.02 
-0.51 ± 0.01 
-0.69 ± 0.03 

-0.901 ± 0.003 
-0.926 ± 0.004 
-0.98 ± 0.04 
-0.869 ± 0.003 

Asian 

Price 

Delta 

Gamma 

Vega 

-0.50 ± 0.02 
-0.49 ± 0.03 
-0.53 ± 0.05 
-0.51 ± 0.02 

-0.70 ± 0.03 
-0.59 ± 0.02 
-0.49 ± 0.03 
-0.64 ± 0.04 

-0.85 ± 0.01 
-0.61 ± 0.03 
-0.50 ± 0.03 
-0.75 ± 0.01 

Double KO 

Price 

Delta 

Gamma 

Vega 

-0.49 ± 0.03 
-0.49 ± 0.02 
-0.45 ± 0.03 
-0.50 ± 0.03 

-0.49 ± 0.02 
-0.52 ± 0.03 
-0.51 ± 0.02 
-0.53 ± 0.02 

-0.56 ± 0.03 
-0.55 ± 0.02 
-0.61 ± 0.02 
-0.57 ± 0.04 

Cliquet 

Price 

Vega 

-0.51 ± 0.02 
-0.48 ± 0.02 

-1.00 ± 0.03 
-0.86 ± 0.04 

-0.72 ± 0.09 
-0.62 ± 0.04 


Table 5: Slopes from linear regression with their errors, as in Table |U In a few cases we show three 
decimals since the MC error is lower. 




(a) Price (b) Delta 




(c) Gamma (d) Vega 

Figure 13: European call option price (a) and greeks (b),(c),(d), Log-Log plots of error ejv versus 
number of simulated paths N = 2 P , p = 9,..., 18, D = 32, e = 10 -3 , L = 30 runs: MC+SD with 
antithetic variables (blue), QMC+SD (green), QMC+BBD (magenta). Linear regression lines are also 
shown. 
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(a) Price 


(b) Delta 




(c) Gamma (d) Vega 

Figure 14: Asian call option, e = 5 x 10 -3 . Other details as in Figure [13] 
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(c) Gamma (d) Vega 

Figure 15: Double Knock-out call option. Details as in Figure [Kj] 




(a) Price (b) Vega 

Figure 16: Cliquet option. Details as in Figure fill 
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We observe what follows. 


1. European option (Figure fT3l) : QMC+BBD outperforms other methods, having the highest 
rate of convergence a. and the smallest intercept. QMC+SD for price and vega has higher 
rate of convergence a but also somewhat higher intercepts than MC+SD. It performs 
marginally better for delta and as good as MC+SD for gamma in terms of a values. 

2. Asian option (Figure fl4l) : for price and vega QMC+BBD and QMC+SD have higher a 
than MC+ SD, with QMC+BBD being the most efficient. They also have slightly higher 
a for delta, but both have lower intercepts than MC. For gamma all methods show similar 
convergence. 

3. Double KO option (Figure fT5l) : QMC+BBD has the highest a although its highest value 
a = 0.61 (for gamma) is lower than a for European and Asian options (with an exception 
of gamma of Asian option). Its intercepts for price, delta and gamma also have the lowest 
values among all methods. The QMC+SD is as efficient as MC. 

4. Cliquet option (Figure [TUP : QMC+SD has the highest a , close to 1.0. It also has the 
lowest intercepts among all methods. QMC+BBD has higher a but similar intercepts in 
comparison with MC+SD. 

We stress that slopes and intercepts shown in the previous figures [T3lfT6l do not depend on the 
details of the simulations, in particular the MC seed or the LDS starting point, since we are 
averaging over L = 30 runs. 

In conclusion, QMC+BBD generally outperforms the other methods, except for asian gamma 
where all methods show similar convergence properties and Cliquet option for which QMC+SD 
is the most efficient method. 


4.4 Speed-Up Analysis 

A typical question with Monte Carlo simulation is “how many scenarios are necessary to achieve 
a given precision?”. When comparing two numerical simulation methods, the typical question 
becomes “how many scenarios may I save using method B instead of method A, preserving the 
same precision?”. 

A useful measure of the relative computational performance of two numerical methods is the 
so called speed-up 5* (a) [KMRZ98al. 1PT961J . It is defined as 


c (+') 


(a) 


Nj j \a) 

N?\a) ’ 


(4.6) 


where, in our context, N*\a) is the number of scenarios using the i-th computational method 
(MC+SD, QMC+SD, or QMC+BBD) needed to reach and maintain a given accuracy a w.r.t. ex¬ 
act or almost exact results. Thus, the speed-up 5* (a) quantifies the computational gain of 
method i w.r.t. method j. 

The speed-up JV* could be evaluated through direct simulation, but this would be extremely 
computationally expensive. Thus we resort to the much simpler algorithm described in Appendix 

m 

We show in Table [6] the results of the speed-up analysis obtained for all methods and for all 
option types described in the previous sections. The speed-up measure clearly shows the relative 
efficiencies of the methods considered for each case. In general, QMC+BBD largely outperforms 
the other methods, with a speed-up factor up to 10 3 (European and Barrier gamma) and a few 
exceptions (Asian delta and gamma, Cliquet). QMC+SD is the best method for Cliquet. We 
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Payoff 

Function 

QMC+SD 
vs MC+SD 

a = 1% a = 0.1% 

QMC+BBD 
vs MC+SD 

a = 1% a = 0.1% 

QMC+BBD 
vs QMC+SD 

a = 1% a = 0.1% 

European 

Price 

3 

6 

30 

140 

10 

20 


Delta 

0.3 

0.5 

20 

100 

60 

200 


Gamma 

0.5 

0.5 

200 

1000 

600 

5000 


Vega 

5 

6 

50 

140 

10 

20 

Asian 

Price 

5 

10 

30 

100 

5 

10 


Delta 

0.2 

0.5 

0.4 

2 

2 

5 


Gamma 

0.5 

- 

0.5 

- 

1 

- 


Vega 

6 

10 

30 

100 

5 

10 

Double KO 

Price 

0.5 

0.8 

5 

10 

10 

15 


Delta 

0.5 

1 

5 

20 

10 

20 


Gamma 

0.7 

1.3 

110 

650 

150 

500 


Vega 

0.5 

0.5 

1.5 

1.5 

3 

3 

Cliquet 

Price 

10 

100 

1 

10 

0.1 

0.1 


Vega 

20 

100 

0.5 

1 

0.02 

0.01 


Table 6: Speed-up S'*(a) of the various numerical methods w.r.t. each other (see columns), for different 
option types. The shift e for finite differences is the same as used in the previous sections. Missing values 
of S* mean that the required accuracy cannot be reached since it is smaller than the bias. 


notice in particular that, in most cases, a ten-fold increase of the accuracy a results in a two-fold 
increase of speed-up 5*(a). However, in a few cases (gamma for European and Cliquet options), 
such an increase can result in up to ten-folds increase of 5*(a). 

The difficulty with speed-up is the possible non-monotonicity of the convergence plot for a 
given numerical methods. Unfortunately, our algorithm to estimate speed-up in Appendix [Bl 
cannot capture unexpected fluctuations of the convergence plot, which could lead to underesti¬ 
mate iV*(a). However, we believe that the choice of the 3-sigma confidence interval in eq. (IB. II) 
makes our speed-up analysis reliable, at least when coupled with the stability analysis described 
in the next Section S3J 

4.5 Stability Analysis 

We have already observed that QMC convergence is often smoother than MC (see Figures 151 1121) : 
such nronotonicity and stability guarantee better convergence for a given number of paths N. In 
order to quantify monotonicity and stability of the various numerical techniques, the following 
strategy is used: we divide the range of path simulations N in 10 windows of equal length, 
and we compute the sample mean rm and the sample standard deviation (“volatility”) s, for 
each window i. Then, log-returns log(mi/mj_i) and volatilities Sj, for i = 2,..., 10, are used 
as measures of, respectively, monotonicity and stability: “monotonic” convergence will show 
non oscillating log-returns converging to zero, “stable” convergence will show low and almost 
flat volatility. We performed stability analysis for MC and QMC methods. For QMC we used 
two different generators: pure QMC with Broda generator and randomized Quasi Monte Carlo 
(rQMC) with Matlab generator The results are shown in Figures [T71I201 


8 Matlab Function sobolset with the MatousekAff ineOwen scrambling method was used. 
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(b) Delta 





(c) Gamma 



Figure 17: Log-returns (upper plots) and volatilities (lower plots) of European call option price (a) and 
greeks (b),(c),(d), for D = 32, e = 10~ 3 , MC+SD (blue), rQMC+BBD (green) and pure QMC+BBD 
(red). The number of simulation paths ranges from 100 to 10,000 grouped in 10 windows each containing 
10 samples (x-axis). 
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(c) Gamma (d) Vega 

Figure 19: Double Knock-out call option. Details as in Figure [lH] 


xIO -3 xUf 3 






(a) Price (b) Vega 

Figure 20: Cliquet option. QMC and rQMC with SD were used here. Other details as in Figure [18l 
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(a) Price 


(b) Delta 





(c) Gamma 



Figure 21: Asian call option with D = 252, e = 5 x 10 3 . Results are shown for rQMC+SD+Matlab 
(green) and rQMC+BBD+Matlab (magenta), and QMC+SD+Broda (blue) and QMC+BBD+Broda 
(red). 
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We observe that, in general, QMC+Broda and rQMC+Matlab are more monotonic and 
stable than MC+SD. However, this fact is less evident for Asian delta and gamma, where QMC 
lacks monotonicity and stability w.r.t. MC, with QMC+BRODA being slightly more stable than 
rQMC+Matlab. As we know from the results of GSA for this case, higher order interactions are 
present and the effective dimensions are large (see Table [3]). 

In order to understand also the effect of dimension D on monotonicity and stability, we 
run a similar experiment for an Asian option with D = 252 fixing dates using both QMC and 
rQMC with SD and BBD. The results are shown in Figure [21J We observe that pure QMC with 
Broda generator preserves monotonicity and stability much more than randomized QMC based 
on Mat lab generator for all cases including delta and gamma, with QMC+BBD+Broda showing 
the best stability. It is also interesting to note that the increase in dimension resulted in the 
decrease in the effective dimensions for the case of the BBD (but not for the SD). 

We conclude that good high-dimensional LDS generators are crucial to obtain a smooth 
monotonic and stable convergence of the Monte Carlo Simulation in high effective dimensional 
problems. 

5 Conclusions 

In this work we presented an updated overview of the application of Quasi Monte Carlo (QMC) 
and Global Sensitivity Analysis (GSA) methods in Finance, w.r.t. standard Monte Carlo (MC) 
methods. In particular, we considered prices and greeks (delta, gamma, vega) for selected 
payoffs with increasing degree of complexity and path-dependency (European Call, Geometric 
Asian Call, Double Barrier Knock-Out, Cliquet options). We compared standard discretization 
(SD) vs Brownian bridge discretization (BBD) schemes of the underlying stochastic diffusion 
process, and different sampling of the underlying distribution using pseudo random vs high 
dimensional Sobol’ low discrepancy sequences. We applied GSA and we performed detailed 
and systematic analysis of convergence diagrams, error estimation, performance, speed-up and 
stability of the different MC and QMC simulations. 

The GSA results in Section 1+21 revealed that effective dimensions associated to QMC+BBD 
simulations are generally lower than those associated to MC+SD simulations, and how much 
such dimension reduction acts for different payoffs and greeks (Figures |T]|S] and Tables DlCTli . 
Effective dimensions, being linked with the structure of ANOVA decompositions (the number of 
important inputs, importance of high order interactions) fully explained the superior efficiency 
of QMC+BBD due to the specifics of Sobol’ sequences and BBD. The BBD is generally more 
efficient than SD, but with some exceptions, Cliquet options in particular. 

The performance analysis results in Section PPl showed that QMC+BBD outperforms MC+SD 
in most cases, showing faster and more stable convergence to exact or almost exact results (Fig¬ 
ures I9KT21 [T3lfT6l and Tables 141151) . with some exceptions such as Asian option gamma where all 
methods showed similar convergence properties. 

The speed-up analysis results in Section 14.41 confirmed that the superior performance of 
QMC+BBD allows significative reduction of the number of scenarios to achieve a given accuracy, 
leading to significative reduction of computational effort (Table [H]) • The size of the reduction 
scales up to 10 3 (European and Double KO gamma), with a few exceptions (Asian delta and 
gamma, Cliquet). 

Finally, the stability analysis results in 14.51 confirmed that QMC+BBD simulations are gen¬ 
erally more stable and monotonic than MC+SD, with the exception of Asian delta and gamma 

(Figures H3EH) • 

We conclude that the methodology presented in this paper, based on Quasi Monte Carlo, 
high dimensional Sobol’ low discrepancy generators, efficient discretization schemes, global sen¬ 
sitivity analysis, detailed convergence diagrams, error estimation, performance, speed-up and 
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stability analysis, is a very promising technique for more complex problems in finance, in partic¬ 
ular, credit/debt/funding/capital valuation adjustments (CVA/DVA/FVA/KVA) and market 
and counterparty risk measured, based on multi-dimensional, multi-step Monte Carlo simula¬ 
tions of large portfolios of trades. Such simulations can run, in typical real cases, ~ 10 2 time 
simulation steps, ~ 10 3 (possibly correlated) risk factors, ~ 10 3 — 10 4 MC scenarios, ~ 10 4 — 10 5 
trades, 60 years maturity, leading to a nominal dimensionality of the order D 10 5 , and to a 
total of 10 9 * — 10 11 evaluations. Unfortunately, a fraction ~ 1% of exotic trades may require dis¬ 
tinct MC simulations for their evaluation, nesting another set of 10 3 — 10 5 MC scenarios, thus 
leading up to 10 14 evaluations. Finally, hedging CVA/DVA/FVA/KVA valuation adjustments 
w.r.t. to their underlying risk factors (typically credit/funding curves) also requires the com¬ 
putation of their corresponding greeks w.r.t. each term structure node, adding another ~ 10 2 
simulations. This is the reason why the industry is continuously looking for advanced tech¬ 
niques to reduce computational times: grid computing, GPU computing, adjoint algorithmic 
differentiation (AAD), etc. (see e.g. |Shel5] h 

We argue that, using QMC sampling (instead of MC) to generate the scenarios of the un¬ 
derlying risk factors and to price exotic trades may significantly improve the accuracy, the 
performance and the stability of such monster-simulations, as shown by preliminary results on 
real portfolios }BKSl~4j . Furthermore, GSA should suggest how to order the risk factors accord¬ 
ing to their relative importance, thus reducing the effective dimensionality. Such applications 
will need further research. 

Appendices 

A Error Optimization in Finite Difference Approximation 

There are two contributions to the root mean square error when greeks are computed by 
MC/QMC simulation via finite differences: variance and bias )Gla03| . The first source of un¬ 
certainty comes from the fact that we are computing prices through simulation over a finite 
number of scenarios, while the latter is due to the approximation of a derivative with a finite 
difference. In order to minimize the variance, we use the same set of (quasi)random numbers for 
the computation of V(6), V(9 + h) and V(9 — h ), where V is the option price, the parameter 6 
is the spot for delta and gamma or the volatility for vega and h is the increment on 0. In order 
to minimize the bias of the finite differences we use central differences, so that it is of the order 
h?. The increments h are chosen to be h = eSo, for A and T, and h = e, for V, for a given “shift 
parameter” e. The choice of the appropriate e is guided by the following considerations. The 
MC/QMC root mean square error estimate of finite differences is given by |Gla03] : 

£ = V / ^ +62 '‘ 4 "' (a ' 1) 

The first term in the square root is a “statistical” error related to the variance c. It depends on 
N as well as on e. a = 0.5 for MC and, usually, 0.5 < a < 1 for QMC, while (3 = 1 for first 
derivatives and (3 = 3 for second derivatives. The second term is the systematic error due to the 
bias of finite differences: it is independent of N but it depends on e. The constant b is given 
by b = g^r(^) f° r central differences of the first order (delta and vega) and b = j^f^r($) for 
central differences of the second order (gamma). One can see that, when h is decreasing, the 
bias term decreases as well while the variance term increases, therefore we fine tune h in such 

9 Some of these metrics, such as EPE/ENE or expected shortfall, are defined as means or conditional means, 

while some other metrics, such as VaR or PFE, are defined as quantiles of appropriate distributions. 
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a way that the variance term is not too high in the relevant range for IV, while the bias term 
remains negligible so that (1 A. 1 1) follows approximately a power law. We note that the optimal 
value of h is not observed to vary too much with N in the range used for our tests. Indeed, it 
can be computed analytically from (1A.1D as the minimum of e: 


i ( Pc \ e+* 
lN \4b 2 N 2a ) 


(A.2) 


We see that the powers l and ^ (corresponding to /3 = 1 and (3 = 3 respectively) largely flatten 
hj\f as a function of N. 


B Speed-Up Computation 

We identify the number of scenarios N*\a) in eq. (14.61) using the i-th computational method 
needed to reach and maintain a given accuracy a as the first number of simulated paths such 
that, for any N > IV* 

V — a < Vn ± 3e < V + a , (B.l) 

where V and Vn are respectively the exact and simulated values of prices or greeks and s is the 
standard error. The threshold IV* could be evaluated through direct simulation, but this would 
be extremely computationally expensive. Extracting IV* from plots defined by (12.23[) can’t be 
applied directly because, in the case of greeks, such plots are correct only for a limited range of 
values of IV, i.e. as long as the bias term in (IA.1D does not become dominant. Extrapolating IV* 
from plots to high values of N is necessary to compute speed-up, but the relation between RMSE 
and N would not be linear anymore. We therefore follow a different procedure to determine IV*. 
Equation (IA.1I) can be rewritten as 


log e = k — a log IV, (B.2) 

where k = ^ log and a are, respectively, the intercept and the slope computed from linear 
regressions on £n given by (12.241 ) . Therefore, IV* is found by imposing 


a = 3 



+ b 2 h 4 , 


(B.3) 


and is given by 


IV* (h, a) 


9 e 2kh \ 
a 2 - 9 b 2 h A ) 



(B.4) 


We have written kh and ah in order to stress that they also depend on the choice of h made 
while carrying out the tests in Section 2) this dependence on h can be stronger than the ex¬ 
plicit dependence in (IB.41) . Constant b is computed from derivatives of V (see discussion after 
presenting equation 1 A. II) ); k and a are the intercepts and slopes correspondingly taken from 
plots (Figures fUlfUD : this is possible, since these plots are obtained for a range of N such that 
the second term in 1 A. II) is negligible. It is clear that the domain of IV* is limited to a > 3 bh 2 . 
In the case of prices, equation (IB,4|) simplifies to 


IV* (a) = 



(B.5) 
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